In this analysis different portfolios consisting of stocks that are among the top 10 constituents of the MSCI Spain Index are backtested using the functionality of the portvine package. First the most important packages are loaded and after that the data will be imported and discussed shortly.

library(portvine)
library(tidyverse)
library(patchwork)
# utility color vector for visualizations
yes_no_cols <- c("#92B8DE", "#db4f59")
theme_set(
  theme_minimal() +
  theme(plot.title = ggtext::element_markdown(size = 11),
        plot.subtitle = ggtext::element_markdown(size = 9))
)
# load the data
load(here::here("data", "msci_spain_data_clean.RData"))

A first glimpse at the data

glimpse(msci_spain_data)
## Rows: 3,089
## Columns: 11
## $ date             <dttm> 2010-01-05, 2010-01-06, 2010-01-07, 2010-01-08, 2010~
## $ msci_spain_index <dbl> 7.227127e-03, -3.519919e-03, -3.271415e-03, -7.783202~
## $ iberdrola        <dbl> 0.0011983863, -0.0013316290, -0.0064166196, -0.001046~
## $ banco_santander  <dbl> 0.010994558, 0.007111767, -0.004592700, 0.003339731, ~
## $ inditex          <dbl> -0.0066202730, 0.0004579803, -0.0147604156, 0.0034790~
## $ cellnex_telecom  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N~
## $ repsol_ypf       <dbl> 0.0029173273, -0.0045062036, 0.0005267713, -0.0010632~
## $ ferrovial        <dbl> 0.013137950, 0.032525280, -0.005260778, 0.003207450, ~
## $ amadeus_it_group <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N~
## $ telefonica       <dbl> -0.0015128243, -0.0073567397, -0.0099731076, -0.01763~
## $ bbv_argentaria   <dbl> 0.0069158977, 0.0038205007, -0.0045843707, 0.00763226~

One can see that there are 11 columns. The date column gives obviously the date as daily return data will be analyzed here. The daily log returns of the overall index are given in the column msci_spain_index and all other columns contain the daily log returns of the respective stock. One can also detect missing values which can be nicely observed in the visualization below.

This makes total sense as for example Cellnex Telecom stocks are traded since 2015, so there can not be earlier return data. So one could now either impute the values in order to use the whole time frame but here the more simplistic approach of just looking at a time frame without any missing values is applied. This reduced data set is available as msci_spain_complete_data and contains only complete cases.

summary(msci_spain_complete_data$date)
##                  Min.               1st Qu.                Median 
## "2015-05-08 00:00:00" "2016-12-21 18:00:00" "2018-08-07 12:00:00" 
##                  Mean               3rd Qu.                  Max. 
## "2018-08-07 16:52:04" "2020-03-24 06:00:00" "2021-11-08 00:00:00"
nrow(msci_spain_complete_data)
## [1] 1696

This means that the time frame is from the 8th Mai of 2015 until the 8th November of 2021 and there are 1696 observations. Below one can have a look at daily log returns of the overall MSCI Spain index.

So one can observe greater volatility during the stock market selloffs 2015-2016 maybe due to Chinese stock market turbulences, the EU dept crisis and the Brexit votum as well as for the first ‘Covid-19 year’ 2020. This should give the opportunity to access how robust the risk measure estimates are during a crisis. With a training data window size of 1000 one would start estimating risk measures from March 2019 so the high volatility corona situation will be covered.

# first risk estimate date for 1000 as the training window size
msci_spain_complete_data$date[1001]
## [1] "2019-03-08 UTC"

Now specify the first portfolio and analyze it unconditionally as well as conditionally.

The missing leader portfolio missl

Here one picks as the stock portfolio equally weighted all the stocks without the one with the highest market capitalization that makes roughly 17% of the MSCI Spain i.e. Iberdrola, that is mainly in the energy business. First one models this portfolio and estimates the risk estimates unconditionally and then conditionally on the biggest player Iberdrola.

Unconditional missl portfolio

First specify suitable marginal model settings as well as settings for the vine copula.

missl_marginal_settings <- marginal_settings(
  train_size = 1000,
  refit_size = 50
)

missl_vine_settings <- vine_settings(
  train_size = 200,
  refit_size = 25,
  family_set = "all",
  vine_type = "rvine"
)

To compare run times one estimates the risk once sequentially and once in parallel.

# sequentially
missl_risk_roll_seq <- estimate_risk_roll(
  msci_spain_complete_data %>% select(-c(date, msci_spain_index, iberdrola)),
  weights = NULL,
  marginal_settings = missl_marginal_settings,
  vine_settings = missl_vine_settings,
  alpha = c(0.01, 0.05, 0.1),
  risk_measures = c("VaR", "ES_mean", "ES_median", "ES_mc"),
  n_samples = 1000,
  n_mc_samples = 1000,
  trace = TRUE
)
## The last window of interest is shorter (width:  46 ) than the specified window width of 50
## 
## Fit marginal models:
## banco_santander  inditex  cellnex_telecom  repsol_ypf  ferrovial  amadeus_it_group  telefonica  bbv_argentaria  
## 
## Fit vine copula models and estimate risk.
## Vine windows:
## (1/28) (2/28) (3/28) (4/28) (5/28) (6/28) (7/28) (8/28) (9/28) (10/28) (11/28) (12/28) (13/28) (14/28) (15/28) (16/28) (17/28) (18/28) (19/28) (20/28) (21/28) (22/28) (23/28) (24/28) (25/28) (26/28) (27/28) (28/28)
# 4 parallel R session
future::plan("multisession", workers = 4)
missl_risk_roll <- estimate_risk_roll(
  msci_spain_complete_data %>% select(-c(date, msci_spain_index, iberdrola)),
  weights = NULL,
  marginal_settings = missl_marginal_settings,
  vine_settings = missl_vine_settings,
  alpha = c(0.01, 0.05, 0.1),
  risk_measures = c("VaR", "ES_mean", "ES_median", "ES_mc"),
  n_samples = 1000,
  n_mc_samples = 1000,
  trace = TRUE
)
## The last window of interest is shorter (width:  46 ) than the specified window width of 50
## 
## Fit marginal models:
## banco_santander  inditex  cellnex_telecom  repsol_ypf  ferrovial  amadeus_it_group  telefonica  bbv_argentaria  
## 
## Fit vine copula models and estimate risk.
## Vine windows:
## (1/28) (2/28) (3/28) (4/28) (5/28) (6/28) (7/28) (8/28) (9/28) (10/28) (11/28) (12/28) (13/28) (14/28) (15/28) (16/28) (17/28) (18/28) (19/28) (20/28) (21/28) (22/28) (23/28) (24/28) (25/28) (26/28) (27/28) (28/28)
future::plan("sequential")

The results and their respective run times can be seen below.

missl_risk_roll_seq
## An object of class <portvine_roll>
## Number of ARMA-GARCH/ marginal windows: 14 
## Number of vine windows: 28 
## Risk measures estimated: VaR ES_mean ES_median ES_mc 
## Alpha levels used: 0.01 0.05 0.1 
## 
## Time taken: 6.5192 minutes
missl_risk_roll
## An object of class <portvine_roll>
## Number of ARMA-GARCH/ marginal windows: 14 
## Number of vine windows: 28 
## Risk measures estimated: VaR ES_mean ES_median ES_mc 
## Alpha levels used: 0.01 0.05 0.1 
## 
## Time taken: 2.1944 minutes
# thus the speedup is already for 4 workers
missl_risk_roll_seq@time_taken / missl_risk_roll@time_taken
## [1] 2.970871
# besides the print method one can also call the summary for more details
summary(missl_risk_roll)
## An object of class <portvine_roll>
## 
## --- Marginal models ---
## Number of ARMA-GARCH/ marginal windows: 14 
## Train size:  1000 
## Refit size:  50 
## 
## --- Vine copula models ---
## Number of vine windows: 28 
## Train size:  200 
## Refit size:  25 
## Vine copula type:  rvine 
## Vine family set:  all 
## 
## --- Risk estimation ---
## Risk measures estimated: VaR ES_mean ES_median ES_mc 
## Alpha levels used: 0.01 0.05 0.1 
## Number of estimated risk measures: 8352 
## Number of samples for each risk estimation: 1000 
## 
## Time taken: 2.1944 minutes.

Inspect marginal and vine copula model fits

Now one can use the accessor functions to extract the fitted marginal as well as the vine copula models. But first one has a close look at the marginal models and can assess their model quality.

### get the fitted marginal models
missl_fitted_marginals <- fitted_marginals(missl_risk_roll)
# for each asset one now can access the uGARCHroll object that for example 
# contains the fitted coefficients for the second rolling window 

# For now analyze the residuals (mean equation) and the squared residuals 
# (volatilitty equation) of every first of each asset.
missl_resid_plots <- sapply(
  setdiff(names(msci_spain_complete_data),
          c("msci_spain_index", "iberdrola", "date")),
  function(asset_name) {
    # use again a utility function from the portvine package
    first_model_resid <- roll_residuals(
      missl_fitted_marginals[[asset_name]], 1
    )
    
    simple_exploratory <- tibble(resid = first_model_resid) %>%
      tibble::rowid_to_column(var = "id") %>%
      ggplot(aes(x = id, y = resid)) +
      geom_line(size = 0.2) +
      labs(x = "t", y = expression(z[t]),
           title = str_to_title(str_replace_all(asset_name, "_", " "))) +
      theme(axis.text.x = element_blank(),
            axis.ticks.x = element_blank())
    
    acf_plot <- tibble(
      acf = as.numeric(acf(first_model_resid, type = "cor", lag.max = 20,
                           plot = FALSE)$acf),
      lag = 0:20
    ) %>%
      filter(lag != 0 & lag <= 10) %>%
      ggplot() +
      geom_hline(yintercept = 0, col = "black", size = 0.3) +
      geom_hline(yintercept = qnorm(c(0.025, 0.975)) /
                   sqrt(length(first_model_resid)),
                 linetype = "longdash", col = yes_no_cols[1], size = 0.5) +
      geom_segment(aes(x = lag, xend = lag, y = 0, yend = acf)) +
      geom_point(aes(x = lag, y = acf)) +
      scale_x_continuous(breaks = seq(1, 10, 1)) +
      ylim(-1, 1) +
      labs(x = "h", y = "ACF(h)")
    
    ljungplot <- tibble(
      pval = sapply(
        1:10,
        function(i) Box.test(first_model_resid, lag = i,
                             type = "Lju")$p.value),
      lag = 1:10) %>%
      ggplot() +
      geom_hline(yintercept = 0, col = "black", size = 0.3) +
      geom_hline(yintercept = 0.05,
                 linetype = "longdash", col = yes_no_cols[1], size = 0.5) +
      geom_line(aes(x = lag, y = pval)) +
      geom_point(aes(x = lag, y = pval)) +
      scale_x_continuous(breaks = seq(1, 10, 1)) +
      labs(x = "h", y = "p-value of Ljung-Box test at lag h")
    
    (simple_exploratory / (ljungplot + acf_plot)) +
      plot_layout(nrow = 2)
  }, USE.NAMES = TRUE, simplify = FALSE)

missl_squaredresid_plots <- sapply(
  setdiff(names(msci_spain_complete_data),
          c("msci_spain_index", "iberdrola", "date")),
  function(asset_name) {
    # use again a utility function from the portvine package
    first_model_resid <- roll_residuals(
      missl_fitted_marginals[[asset_name]], 1
    )^2
    
    simple_exploratory <- tibble(resid = first_model_resid) %>%
      tibble::rowid_to_column(var = "id") %>%
      ggplot(aes(x = id, y = resid)) +
      geom_line(size = 0.2) +
      labs(x = "t", y = expression(z[t]^2),
           title = str_to_title(str_replace_all(asset_name, "_", " "))) +
      theme(axis.text.x = element_blank(),
            axis.ticks.x = element_blank())
    
    acf_plot <- tibble(
      acf = as.numeric(acf(first_model_resid, type = "cor", lag.max = 20,
                           plot = FALSE)$acf),
      lag = 0:20
    ) %>%
      filter(lag != 0 & lag <= 10) %>%
      ggplot() +
      geom_hline(yintercept = 0, col = "black", size = 0.3) +
      geom_hline(yintercept = qnorm(c(0.025, 0.975)) /
                   sqrt(length(first_model_resid)),
                 linetype = "longdash", col = yes_no_cols[1], size = 0.5) +
      geom_segment(aes(x = lag, xend = lag, y = 0, yend = acf)) +
      geom_point(aes(x = lag, y = acf)) +
      scale_x_continuous(breaks = seq(1, 10, 1)) +
      ylim(-1, 1) +
      labs(x = "h", y = "ACF(h)")
    
    ljungplot <- tibble(
      pval = sapply(
        1:10,
        function(i) Box.test(first_model_resid, lag = i,
                             type = "Lju")$p.value),
      lag = 1:10) %>%
      ggplot() +
      geom_hline(yintercept = 0, col = "black", size = 0.3) +
      geom_hline(yintercept = 0.05,
                 linetype = "longdash", col = yes_no_cols[1], size = 0.5) +
      geom_line(aes(x = lag, y = pval)) +
      geom_point(aes(x = lag, y = pval)) +
      scale_x_continuous(breaks = seq(1, 10, 1)) +
      labs(x = "h", y = "p-value of Ljung-Box test at lag h")
    
    (simple_exploratory / (ljungplot + acf_plot)) +
      plot_layout(nrow = 2)
  }, USE.NAMES = TRUE, simplify = FALSE)
# exemplary 2 plots, all plots showed basically the same behavior like the ones
# below
missl_resid_plots$telefonica

missl_squaredresid_plots$telefonica

missl_resid_plots$repsol_ypf

missl_squaredresid_plots$repsol_ypf

# utility function to visualize the pvalues of the Ljung box tests for all
# assets simultaneously via a heatmap
ljung_heatmap <- function(roll, roll_num = 1) {
  asset_names <- fitted_vines(roll)[[1]]$names
  roll_marginals <- fitted_marginals(roll)
  ljung_data <-sapply(asset_names, function(asset_name) {
    # use again a utility function from the portvine package
    model_resid <- roll_residuals(
      roll_marginals[[asset_name]], roll_num = roll_num
    )
    sapply(
        1:10,
        function(i) Box.test(model_resid, lag = i, type = "Lju")$p.value
    )
  }, USE.NAMES = TRUE, simplify = TRUE)
  ljung_data %>%
    as_tibble() %>%
    rowid_to_column("lag") %>%
    pivot_longer(-lag, names_to = "asset", values_to = "pval") %>%
    ggplot(aes(x = lag, y = asset, fill = pval)) +
    geom_tile() +
    scale_y_discrete(labels = ~ str_to_title(str_replace_all(.x, "_", " "))) +
    scale_x_continuous(breaks = 1:10) +
    labs(y = "", x = "h", fill = "p-value",
         title = "Results of the Ljung-Box tests",
         caption = paste("Rolling window:", roll_num)) +
    scale_fill_gradientn(colours = c(yes_no_cols[2],"#C37285" ,
                                     yes_no_cols[1], "#2a82db"),
                           values = scales::rescale(c(0, 0.05 - 0.01, 
                                                          0.05, 1)), 
                           breaks = c(0.05), 
                           labels = c(0.05),
                           guide = guide_colourbar(nbin = 1000)) +
    theme(legend.position = "right",
          panel.grid.minor.x = element_blank())
}
# use the function:
ljung_heatmap(missl_risk_roll)

ljung_heatmap(missl_risk_roll, 3)

So overall the marginal models seem to satisfy the model assumptions. Notably there were many cases with a heavier lower tail which can be explained by the high volatility corona crisis. As already the skewed student’s t distribution was chosen for the residual distribution one already accounts for heavy tails.

Now one can have a look at the fitted vine copulas. Again there is an appropriate accessor function that helps with extracting them. After that one can visually inspect how the first tree of the vines has changed over time.

### get the fitted vine copula models
missl_fitted_vines <- fitted_vines(missl_risk_roll)
# then one can for example have a look at how the fitted vine structure evolved
#  over time by plotting the first tree for different rolling windows

# utility function for labeled vinecop tree plots
labeled_vinecop_plot <- function(vine) {
  node_labels <- vine$names[vine$structure$order]
  ff <- tempfile()
  png(filename = ff)
  plot <- rvinecopulib:::plot.vinecop(vine) 
  dev.off()
  unlink(ff)
  plot$data$name <- str_to_title(str_replace_all(node_labels, "_", " "))
  plot
    # ggraph::geom_node_label(
    #   aes(label = str_to_title(str_replace_all(node_labels, "_", " "))) 
    # )
}
# utility for table of used bivariate copula families
bicops_used <- function(vine) {
  table(unlist(vine$pair_copulas)[names(unlist(vine$pair_copulas)) == "family"])
}

labeled_vinecop_plot(missl_fitted_vines[[1]])

bicops_used(missl_fitted_vines[[1]])
## 
##      bb1      bb8  clayton    frank gaussian   gumbel    indep      joe 
##        1        1        2        8        4        5        5        2
labeled_vinecop_plot(missl_fitted_vines[[14]])

bicops_used(missl_fitted_vines[[14]])
## 
##  clayton    frank gaussian   gumbel    indep      joe        t 
##        3        2        4        8        7        2        2
labeled_vinecop_plot(missl_fitted_vines[[28]])

bicops_used(missl_fitted_vines[[28]])
## 
##      bb1      bb7  clayton    frank gaussian   gumbel    indep      joe 
##        1        1        1        8        3        2        7        3 
##        t 
##        2

The first two assets seem to be very important in their dependence structure.

Analyze estimated risk measures

Now finally the risk estimates are going to be analyzed. First one compares the 3 different methods of estimating the Expected Shortfall.

Also one might compare the number of exceedances.

risk_estimates(missl_risk_roll,
               risk_measures = c("ES_mean", "ES_median", "ES_mc"),
               alpha = 0.05,
               exceeded = TRUE) %>%
  group_by(risk_measure) %>%
  summarise(relative_exceedances = mean(exceeded))

Here actually we find an indication for the fact that the values that fell below the corresponding VaR where left skewed which leads to the fact that the mean will be a more conservative estimate in those cases which might explain the slightly lower exceedance rate.

Exceedance plots

risk_estimates(missl_risk_roll,
               risk_measures = c("ES_mean", "VaR"),
               alpha = 0.01) %>%
  ggplot() +
  geom_line(aes(x = row_num, y = realized), col = "lightgrey") +
  geom_line(aes(x = row_num, y = risk_est,
                col = risk_measure, linetype = risk_measure),
            size = 0.5) +
  labs(x = "estimation window", y = "portfolio log returns",
       linetype = "Risk measure",
       title = "Comparison of risk measure behaivior for alpha level 0.01") +
  scale_color_manual(values = c(yes_no_cols[1], "#477042"),
                     name = "Risk measure") +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

risk_estimates(missl_risk_roll,
               risk_measures = c("VaR"),
               alpha = 0.01,
               exceeded = TRUE) %>%
  ggplot() +
  geom_line(aes(x = row_num, y = realized), col = "lightgrey") +
  geom_line(aes(x = row_num, y = risk_est), col = yes_no_cols[1]) +
  geom_point(aes(x = row_num, y = realized), col = yes_no_cols[2], 
             inherit.aes = FALSE, data = . %>% filter(exceeded)) +
  labs(x = "estimation window",
       y = "portfolio log returns",
       col = "Exceeded",
       subtitle = paste0("Exceedances are highlighted in ",
                     "<span style='color:",
                     yes_no_cols[2],
                     "'>**red**</span>",
                     "."),
       title = "Risk measure: VaR with alpha level 0.01"
       ) +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "none")

risk_estimates(missl_risk_roll,
               risk_measures = c("ES_mean"),
               alpha = 0.01,
               exceeded = TRUE) %>%
  ggplot() +
  geom_line(aes(x = row_num, y = realized), col = "lightgrey") +
  geom_line(aes(x = row_num, y = risk_est), col = yes_no_cols[1]) +
  geom_point(aes(x = row_num, y = realized), col = yes_no_cols[2], 
             inherit.aes = FALSE, data = . %>% filter(exceeded)) +
  labs(x = "estimation window",
       y = "portfolio log returns",
       col = "Exceeded",
       subtitle = paste0("Exceedances are highlighted in ",
                     "<span style='color:",
                     yes_no_cols[2],
                     "'>**red**</span>",
                     "."),
       title = "Risk measure: ES (mean estimation) with alpha level 0.01"
       ) +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "none")

Also fit the unconditional missl portfolio with only D vine copulas allowed to compare these D vines with the then conditional ones.

missl_vine_settings_d <- vine_settings(
  train_size = 200,
  refit_size = 25,
  family_set = "all",
  vine_type = "dvine"
)
future::plan("multisession", workers = 6)
missl_risk_roll_d <- estimate_risk_roll(
  msci_spain_complete_data %>% select(-c(date, msci_spain_index, iberdrola)),
  weights = NULL,
  marginal_settings = missl_marginal_settings,
  vine_settings = missl_vine_settings_d,
  alpha = c(0.01, 0.05, 0.1),
  risk_measures = c("VaR", "ES_mean", "ES_median", "ES_mc"),
  n_samples = 1000,
  n_mc_samples = 1000,
  trace = TRUE
)
## The last window of interest is shorter (width:  46 ) than the specified window width of 50
## 
## Fit marginal models:
## banco_santander  inditex  cellnex_telecom  repsol_ypf  ferrovial  amadeus_it_group  telefonica  bbv_argentaria  
## 
## Fit vine copula models and estimate risk.
## Vine windows:
## (1/28) (2/28) (3/28) (4/28) (5/28) (6/28) (7/28) (8/28) (9/28) (10/28) (11/28) (12/28) (13/28) (14/28) (15/28) (16/28) (17/28) (18/28) (19/28) (20/28) (21/28) (22/28) (23/28) (24/28) (25/28) (26/28) (27/28) (28/28)
future::plan("sequential")

Conditional missl portfolio

Estimation is performed in the same fashion as above with some additional arguments concerning the conditioning on the stock Iberdrola.

missl_cond_marginal_settings <- marginal_settings(
  train_size = 1000,
  refit_size = 50
)

missl_cond_vine_settings <- vine_settings(
  train_size = 200,
  refit_size = 25,
  family_set = "parametric",
  vine_type = "dvine"
)

future::plan("multisession", workers = 6)
missl_cond_risk_roll <- estimate_risk_roll(
  msci_spain_complete_data %>% select(-c(date, msci_spain_index)),
  weights = NULL,
  marginal_settings = missl_cond_marginal_settings,
  vine_settings = missl_cond_vine_settings,
  alpha = c(0.01, 0.05, 0.1),
  risk_measures = c("VaR", "ES_mean", "ES_median", "ES_mc"),
  n_samples = 1000,
  n_mc_samples = 1000,
  cond_vars = "iberdrola",
  cond_alpha = c(0.05, 0.5),
  trace = TRUE
)
## The last window of interest is shorter (width:  46 ) than the specified window width of 50
## 
## Fit marginal models:
## iberdrola  banco_santander  inditex  cellnex_telecom  repsol_ypf  ferrovial  amadeus_it_group  telefonica  bbv_argentaria  
## 
## Fit vine copula models and estimate risk.
## Vine windows:
## (1/28) (2/28) (3/28) (4/28) (5/28) (6/28) (7/28) (8/28) (9/28) (10/28) (11/28) (12/28) (13/28) (14/28) (15/28) (16/28) (17/28) (18/28) (19/28) (20/28) (21/28) (22/28) (23/28) (24/28) (25/28) (26/28) (27/28) (28/28)
future::plan("sequential")
summary(missl_cond_risk_roll)
## An object of class <cond_portvine_roll>
## 
## --- Conditional settings ---
## Conditional variable(s): iberdrola 
## Number of conditional estimated risk measures: 16704 
## Conditioning quantiles: 0.05 0.5 
## 
## --- Marginal models ---
## Number of ARMA-GARCH/ marginal windows: 14 
## Train size:  1000 
## Refit size:  50 
## 
## --- Vine copula models ---
## Number of vine windows: 28 
## Train size:  200 
## Refit size:  25 
## Vine copula type:  dvine 
## Vine family set:  parametric 
## 
## --- Risk estimation ---
## Risk measures estimated: VaR ES_mean ES_median ES_mc 
## Alpha levels used: 0.01 0.05 0.1 
## Number of estimated risk measures: 8352 
## Number of samples for each risk estimation: 1000 
## 
## Time taken: 4.5191 minutes.

A short look at the p-values of the Ljung Box tests.

ljung_heatmap(missl_cond_risk_roll)

Compare the fitted conditional vines with the ones of the unconditional dvine fitting.

# unconditional first
labeled_vinecop_plot(fitted_vines(missl_risk_roll_d)[[1]])

bicops_used(fitted_vines(missl_risk_roll_d)[[1]])
## 
##      bb1      bb8  clayton    frank gaussian   gumbel    indep      joe 
##        1        2        3        7        3        3        5        4
# conditional first
labeled_vinecop_plot(fitted_vines(missl_cond_risk_roll)[[1]])

bicops_used(fitted_vines(missl_cond_risk_roll)[[1]])
## 
##      bb1      bb7      bb8  clayton    frank gaussian   gumbel    indep 
##        1        1        2        1       10        2        5        7 
##      joe        t 
##        6        1
# unconditional last
labeled_vinecop_plot(fitted_vines(missl_risk_roll_d)[[28]])

# conditional last
labeled_vinecop_plot(fitted_vines(missl_cond_risk_roll)[[28]])

# what about the evolution of the conditional values (iberdrola)
risk_estimates(missl_cond_risk_roll, risk_measures = "ES_mean",
               alpha = 0.05, exceeded = TRUE) %>%
  ggplot() +
  geom_line(
    data = msci_spain_complete_data[1001:nrow(msci_spain_complete_data), ],
    aes(x = 1001:nrow(msci_spain_complete_data),
        y = iberdrola), col = "lightgrey", size = .3) +
  geom_line(aes(x = row_num, y = iberdrola, col = factor(cond_alpha)),
            size = 0.5) +
  scale_color_manual(values = c(yes_no_cols[1], "#477042"),
                     name = "Conditional alpha level") +
  labs(title = "Conditional variable: Iberdrola",
       x = "estimation window", y = "log returns") +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

Nice to see so we expect the corresponding conditional risk measures with low conditional \(\alpha\) level to be way more conservative under the assumption that Iberdrola has a generally positive dependence with the portfolio of interest. Therefore look at the next plot.

risk_estimates(missl_cond_risk_roll, risk_measures = "ES_mean",
               alpha = 0.05, exceeded = TRUE) %>%
  ggplot() +
  geom_line(aes(x = row_num, y = realized), col = "lightgrey") +
  geom_line(aes(x = row_num, y = risk_est,
                col = factor(cond_alpha), linetype = factor(cond_alpha)),
            size = 0.5) +
  labs(x = "estimation window", y = "portfolio log returns",
       linetype = "Conditional alpha level",
       title = "Comparison of risk measure behaivior for 2 conditional alpha levels") +
  scale_color_manual(values = c(yes_no_cols[1], "#477042"),
                     name = "Conditional alpha level") +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

A close look wrt exceedences.

risk_estimates(missl_cond_risk_roll, risk_measures = "ES_mean",
               cond_alpha = 0.5,
               alpha = 0.05, exceeded = TRUE) %>%
  ggplot() +
  geom_line(aes(x = row_num, y = realized), col = "lightgrey") +
  geom_line(aes(x = row_num, y = risk_est), col = yes_no_cols[1]) +
  geom_point(aes(x = row_num, y = realized), col = yes_no_cols[2], 
             inherit.aes = FALSE, data = . %>% filter(exceeded)) +
  labs(x = "estimation window",
       y = "portfolio log returns",
       col = "Exceeded",
       subtitle = paste0("Exceedances are highlighted in ",
                     "<span style='color:",
                     yes_no_cols[2],
                     "'>**red**</span>",
                     "."),
       title = "Risk measure: ES (mean estimation) with alpha level 0.05, conditional alpha level 0.5"
       ) +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "none")

risk_estimates(missl_cond_risk_roll, risk_measures = "ES_mean",
               cond_alpha = 0.05,
               alpha = 0.05, exceeded = TRUE) %>%
  ggplot() +
  geom_line(aes(x = row_num, y = realized), col = "lightgrey") +
  geom_line(aes(x = row_num, y = risk_est), col = yes_no_cols[1]) +
  geom_point(aes(x = row_num, y = realized), col = yes_no_cols[2], 
             inherit.aes = FALSE, data = . %>% filter(exceeded)) +
  labs(x = "estimation window",
       y = "portfolio log returns",
       col = "Exceeded",
       subtitle = paste0("Exceedances are highlighted in ",
                     "<span style='color:",
                     yes_no_cols[2],
                     "'>**red**</span>",
                     "."),
       title = "Risk measure: ES (mean estimation) with alpha and conditional alpha level 0.05"
       ) +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "none")

While the 0.5 conditional value leads to lots of exceedances (not surprisingly) the conditional value of 0.05 shows not many exceedances. For a rough comparison one can have a short comparison with the above fitted unconditional ES with mean estimation for the same confidence level.

risk_estimates(missl_risk_roll,
               risk_measures = c("ES_mean"),
               alpha = 0.05,
               exceeded = TRUE) %>%
  ggplot() +
  geom_line(aes(x = row_num, y = realized), col = "lightgrey") +
  geom_line(aes(x = row_num, y = risk_est), col = yes_no_cols[1]) +
  geom_point(aes(x = row_num, y = realized), col = yes_no_cols[2], 
             inherit.aes = FALSE, data = . %>% filter(exceeded)) +
  labs(x = "estimation window",
       y = "portfolio log returns",
       col = "Exceeded",
       subtitle = paste0("Exceedances are highlighted in ",
                     "<span style='color:",
                     yes_no_cols[2],
                     "'>**red**</span>",
                     "."),
       title = "Risk measure: ES (mean estimation) with alpha level 0.05"
       ) +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    legend.position = "none")

So in the direct comparison the conditional risk measure is way more conservative and might help in stresstesting situations for projected hard times of specific conditional assets.

Double conditional short trial

Here just the estimation with two conditional values is tested shortly.

cond2test_marginal_settings <- marginal_settings(
  train_size = 1000,
  refit_size = 50
)

cond2test_vine_settings <- vine_settings(
  train_size = 200,
  refit_size = 25,
  family_set = "parametric",
  vine_type = "dvine"
)

future::plan("multisession", workers = 6)
cond2test_risk_roll <- estimate_risk_roll(
  msci_spain_complete_data %>% select(-c(date, msci_spain_index)),
  weights = NULL,
  marginal_settings = cond2test_marginal_settings,
  vine_settings = cond2test_vine_settings,
  alpha = c(0.01, 0.05, 0.1),
  risk_measures = c("VaR", "ES_mean", "ES_median", "ES_mc"),
  n_samples = 1000,
  n_mc_samples = 1000,
  cond_vars = c("iberdrola", "banco_santander"),
  cond_alpha = c(0.05, 0.5),
  trace = TRUE
)
## The last window of interest is shorter (width:  46 ) than the specified window width of 50
## 
## Fit marginal models:
## iberdrola  banco_santander  inditex  cellnex_telecom  repsol_ypf  ferrovial  amadeus_it_group  telefonica  bbv_argentaria  
## 
## Fit vine copula models and estimate risk.
## Vine windows:
## (1/28) (2/28) (3/28) (4/28) (5/28) (6/28) (7/28) (8/28) (9/28) (10/28) (11/28) (12/28) (13/28) (14/28) (15/28) (16/28) (17/28) (18/28) (19/28) (20/28) (21/28) (22/28) (23/28) (24/28) (25/28) (26/28) (27/28) (28/28)
future::plan("sequential")
summary(cond2test_risk_roll)
## An object of class <cond_portvine_roll>
## 
## --- Conditional settings ---
## Conditional variable(s): iberdrola banco_santander 
## Number of conditional estimated risk measures: 16704 
## Conditioning quantiles: 0.05 0.5 
## 
## --- Marginal models ---
## Number of ARMA-GARCH/ marginal windows: 14 
## Train size:  1000 
## Refit size:  50 
## 
## --- Vine copula models ---
## Number of vine windows: 28 
## Train size:  200 
## Refit size:  25 
## Vine copula type:  dvine 
## Vine family set:  parametric 
## 
## --- Risk estimation ---
## Risk measures estimated: VaR ES_mean ES_median ES_mc 
## Alpha levels used: 0.01 0.05 0.1 
## Number of estimated risk measures: 8352 
## Number of samples for each risk estimation: 1000 
## 
## Time taken: 4.7541 minutes.
risk_estimates(cond2test_risk_roll)